function result = splitPic(filename)
% close all;
% clear;
% filepath = "test2.png";
% [filepath,name,ext] = fileparts(filepath);
% filename = strcat(filepath, "/", name);
I_cam = imread(filename);
%figure;
%imshow(I_cam);
%load TEMPLATEDATA

%I = I_cam(10+(1:460),90+(1:460));
I = I_cam;
picture = rgb2gray(I);
figure;%         U = edge(t,'canny');
%         %[H, T, R] = hough(I,'Theta',-10:0.1:10);
%         [H, T, R] = hough(U);
%         P = houghpeaks(H,1);
%         lines = houghlines(U,T,R,P);
%         for m=1:length(lines)
%             xy = [lines(m).point1;lines(m).point2];
%             t = t.*~xor(t,imfill(t,xy(:,1)));
%         end

subplot(2,2,1);
imshow(picture);
se = strel('rectangle', [3,3]);
Ibot = imbothat(picture, se);
subplot(2,2,2);
imshow(Ibot);
I = Ibot;
%H = fspecial('laplacian', 0);
%hp1 = imfilter(I, H, 'replicate');
%p1 = hp1;
%S = I;
%I = p1;
makebw = @(I) im2bw(I.data,0.2);
I = I < 30;
subplot(2,2,3);
imshow(I);

I = (picture/2 + uint8(I)*127*1.2);


subplot(2,2,4);
imshow(I);
subplot(2,2,3)
makebw = @(I) im2bw(I.data,median(double(I.data(:)))/1.2/255);
I = ~blockproc(I,[92 92],makebw);
%figure;
%imshow(I);
%% Remove Noise
I = bwareaopen(I,20);
%imshow(I);
%% Clear the border
I = imclearborder(I);
%imshow(I);


%% Find the largest box
hold on;
R = regionprops(I,'Area','BoundingBox','PixelList');
NR = numel(R);
%hold off;

maxArea = 0;
for k = 1:NR
    A(k) = prod(R(k).BoundingBox(3:4));
    if R(k).Area > maxArea
        maxArea = R(k).Area;
        kmax = k;
    end
end


BBmax = R(kmax).BoundingBox;
DIAG1 = sum(R(kmax).PixelList,2);
DIAG2 = diff(R(kmax).PixelList,[],2);

[m,dUL] = min(DIAG1);    [m,dDR] = max(DIAG1);
[m,dDL] = min(DIAG2);    [m,dUR] = max(DIAG2);

pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
h_pts = plot(pts(:,1),pts(:,2),'m','linewidth',3);%HUA TU XIANG KUANG

XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];

T = cp2tform(pts(1:4,:),500*[0 0; 1 0; 1 1; 0 1],'projective');
IT = imtransform(double(I),T,'XData',[1 500], 'YData',[1 500]);

I = IT;
I = im2bw(I);
%hold on;
R = regionprops(I,'Area','BoundingBox','PixelList');

%imshow(I)
NR = numel(R);

maxArea = 0;
for k = 1:NR
    A(k) = prod(R(k).BoundingBox(3:4));
    if R(k).Area > maxArea
        maxArea = R(k).Area;
        kmax = k;
    end
end


BBmax = R(kmax).BoundingBox;
DIAG1 = sum(R(kmax).PixelList,2);
DIAG2 = diff(R(kmax).PixelList,[],2);

[m,dUL] = min(DIAG1);    [m,dDR] = max(DIAG1);
[m,dDL] = min(DIAG2);    [m,dUR] = max(DIAG2);

pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
%h_pts = plot(pts(:,1),pts(:,2),'m','linewidth',3);%HUA TU XIANG KUANG

XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];

%% Identify objects inside the box 
A_tmin = 30; % Bounds for the digit pixel area
A_tmax = 3000;
digitbox_minarea = 20; % Bounds for the digit bounding box area
digitbox_maxarea = 50^2;
tsize = 16;
pictures = zeros(2*tsize+1,2*tsize+1,9, 9);
%figure;
%imshow(I);

cutsize = 5;

kgood = zeros(1,NR);
Pnew = zeros(NR,2);
        for k = 1:NR
            if R(k).Area < A_tmax && A(k) > digitbox_minarea && A(k) < digitbox_maxarea ...
                    && R(k).BoundingBox(3) < 40 && R(k).BoundingBox(4) < 40 ...
                    && R(k).BoundingBox(3) > 1 && R(k).BoundingBox(4) > 1
                
                Pnew(k,:) = [R(k).BoundingBox(1)+R(k).BoundingBox(3)/2 R(k).BoundingBox(2)+R(k).BoundingBox(4)/2];
                
                
                if inpolygon(Pnew(k,1),Pnew(k,2),pts(:,1),pts(:,2))
                    %h_digitcircles(k) = plot(Pnew(k,1),Pnew(k,2),'ro','markersize',24);
                    try
                        t = I(fix(Pnew(k,2)-tsize):fix(Pnew(k,2)+tsize),fix(Pnew(k,1)-tsize):fix(Pnew(k,1)+tsize));
                    catch
                        continue
                    end
          
                    xpos = fix(Pnew(k,2)*9 / 500) + 1;
                    ypos = fix(Pnew(k,1)*9 / 500) + 1;
                    if (abs(0.5 - abs((xpos - 1 - Pnew(k,2)*9 / 500))) > 0.30 || abs(0.5 - abs(ypos - 1 - Pnew(k,1)*9 / 500)) > 0.3)
                        continue
                    end
                    disp([xpos, ypos]);
                    s = size(t);
                    pictures(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize,xpos, ypos) = t(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize);
%                     figure;
%                     imshow(pictures(:,:,xpos,ypos))
                end
                
            end
        end
        figure;
rsize = 30;

Rawpictures = zeros(2*rsize+1,2*rsize+1,9, 9);
twoBitPics = zeros(28, 28, 9, 9);
zeroPic = imread('zeroPic.png');

for j = 1:9
    for k = 1:9
        subplot(10,9,j*9+k);
        imshow(pictures(:,:,j,k));
        t = I(max(1, fix((j-1)*500/9)+1-10):min(fix((j-1)*500/9+2*rsize-10)+1,500),max(1,fix((k-1)*500/9)+1-10):min(fix((k-1)*500/9+2*rsize)+1-10, 500));
%         U = edge(t,'canny');
%         %[H, T, R] = hough(I,'Theta',-10:0.1:10);
%         [H, T, R] = hough(U);
%         P = houghpeaks(H,1);
%         lines = houghlines(U,T,R,P);
%         for m=1:length(lines)
%             xy = [lines(m).point1;lines(m).point2];
%             t = t.*~xor(t,imfill(t,xy(:,1)));
%         end
        %imshow(t);
        s = size(t);
        Rawpictures(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize,j,k) = t(1+cutsize:s(1)-cutsize,1+cutsize:s(2)-cutsize);
        %imshow(Rawpictures(:,:,j,k));
        
        if sum(pictures(:,:,j,k))==0
            twoBitPics(:,:,j,k) = zeroPic;
        else 
            twoBitPics(:,:,j,k) = imresize(pictures(:,:,j,k), [28 28]);
            thresh = graythresh(twoBitPics(:,:,j,k)); %自动确定阈值
            twoBitPics(:,:,j,k) = im2bw(twoBitPics(:,:,j,k),thresh); %对图像二值化
        end
        
    end
end

result = twoBitPics;

